require(tidyverse)
require(sf)
require(ggthemes)

# load data
load("01_Data/spatial_dimension_data.Rdata")
map <- geojsonsf::geojson_sf("01_Data/landkreise_simplify200.geojson")

#### Figure 1: Main descriptive plot ####
p1 <- ggplot() +
  geom_sf(
    data = map,
    fill = "#EFECE6",
    color = "grey20",
    size = 0.1
  ) +
  geom_sf(
    data = location_level,
    aes(fill = delay),
    color = "grey20",
    shape = 21,
    size = 2
  ) +
  scale_fill_gradientn(
    limits = c(-1100, 1100),
    colors = c(
      "black",
      "darkblue",
      "#47A8BD",
      "white",
      "#F27059",
      "#F25C54" ,
      "darkred"
    )
  ) +
  theme_void() +
  theme(legend.position = "left") +
  labs(fill = "Deviation (Days)") +
  ggtitle("A)") +
  theme(plot.title.position = "plot")+
  theme(text = element_text(size = 15)) 


p12 <- ggplot(model_base, aes(x = delay)) +
  geom_histogram(bins = 60, color = "black", fill = "#EFECE6") +
  xlab("Deviation (Days)") +
  ylab("Count") +
  theme_base()+
  theme(text = element_text(size = 15),
        plot.background = element_blank()) 

plot_df <- cbind(location_level, location_coordinates)

p2 <- ggplot(map) +
  geom_sf(fill = "#EFECE6", size = 0.1, color = "grey20") +
  coord_sf() +
  xlab("") +
  ylab("") +
  theme_void() +
  geom_point(aes(X, Y),
             data = plot_df,
             shape = 21,
             size = 0.8) +
  labs(color = "Distance (Hours)") +
  theme(legend.position = "left") +
  geom_segment(data = plot_df, aes(
    x = X,
    y = Y,
    xend = rp_lon,
    yend = rp_lat,
    color = time
  )) +
  scale_color_viridis_c(direction = -1) +
  geom_point(
    aes(rp_lon, rp_lat),
    data = plot_df,
    fill = "#FFBF69",
    color = "black",
    size = 3,
    alpha = 0.9,
    shape = 23
  ) +
  labs(title = "B)") +
  theme(plot.title.position = "plot")+
  theme(text = element_text(size = 15)) 


p2

p22 <- ggplot(model_base, aes(x = time)) +
  geom_histogram(fill = "#EFECE6", color = "black", bins = 60) +
  xlab("Distance (Hours)") +
  ylab("Count") +
  theme_base()+
  theme(text = element_text(size = 15),
        plot.background = element_blank()) 


p3 <- ggplot() +
  geom_sf(
    data = map,
    fill = "#EFECE6",
    color = "grey20",
    size = 0.1
  ) +
  geom_sf(data = location_level, aes(color = risk)) +
  scale_color_continuous(high = "47A8BD", low = "darkblue",
                         breaks = c(1,2,3),
                         labels = c("high","mediun","low"),
                         trans = "reverse") +
  theme_minimal() +
  theme_void() +
  theme(legend.position = "left") +
  labs(color = "Risk (Level)") +
  ggtitle("C)") +
  theme(plot.title.position = "plot")+
  theme(text = element_text(size = 15))

p3

p33 <- model_base %>% 
  group_by(risk_factor) %>% 
  count() %>% 
  mutate(risk_factor=fct_relevel(risk_factor,"low","medium","high")) %>% 
  ggplot(aes(x = risk_factor,y=n)) +
  geom_col(fill = "#EFECE6", color = "black") +
  xlab("Risk (Level)") +
  ylab("Count") +
  theme_base()+
  theme(text = element_text(size = 15),
        plot.background = element_blank())  

p33

pdf("02_Figures/figure01.pdf", width = 16, height = 20)
gridExtra::grid.arrange(p1,
                        p12,
                        p2,
                        p22,
                        p3,
                        p33,
                        ncol = 2,
                        widths = c(10, 5))
dev.off()

#### Figure A4: Boxplot main independent variables#####

model_base$risk_1<- 0
model_base$risk_1[model_base$risk==1] <- 1

model_base$risk_2<- 0
model_base$risk_2[model_base$risk==2] <- 1

model_base$risk_3<- 0
model_base$risk_3[model_base$risk==3] <- 1

pdf("02_Figures/appendix_A4.pdf", width=10, height=5)
model_base %>% 
  mutate(risk_factor=fct_relevel(risk_factor,"low","medium","high")) %>% 
  ggplot(aes(y=risk_factor, x=time))+
  geom_boxplot(width=.3, fill="grey", alpha=0.2)+
  theme_base()+
  xlab("Distance (Hours)")+
  ylab("Risk (Level)")+
  theme(plot.background = element_blank())
dev.off()




A1 <- paste("r =", round(cor(model_base$time, model_base$risk_1),2))
A2 <- paste("r =", round(cor(model_base$time, model_base$risk_2),2))
A3 <- paste("r =", round(cor(model_base$time, model_base$risk_3),2))

rm(list = ls())

